{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian data analysis\n",
    "##  Chapter 6, demo 2\n",
    "\n",
    "Posterior predictive checking  \n",
    "Binomial example - Testing sequential dependence example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os, sys\n",
    "# add utilities directory to path\n",
    "util_path = os.path.abspath(os.path.join(os.path.pardir, 'utilities_and_data'))\n",
    "if util_path not in sys.path and os.path.exists(util_path):\n",
    "    sys.path.insert(0, util_path)\n",
    "\n",
    "# import from utilities\n",
    "import plot_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# edit default plot settings\n",
    "plt.rc('font', size=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Testing sequential dependence example (Gelman et al p. 163)\n",
    "y = np.array([1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0])\n",
    "Ty = np.count_nonzero(np.diff(y))\n",
    "\n",
    "# sufficient statistics\n",
    "n = len(y)\n",
    "s = y.sum()\n",
    "\n",
    "nsamp = 10000\n",
    "t = np.random.beta(s+1, n-s+1, size=nsamp)\n",
    "yr = np.random.rand(n, nsamp) < t\n",
    "# sadly np.count_nonzero does not (yet) support axis parameter\n",
    "Tyr = (np.diff(yr, axis=0) != 0).sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEUCAYAAABkhkJAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuYHFWd//H3hwkhwAQSiISgYCIoiUASIGGjrCYoIEQT\nFlGWBUxE/S3gFQExoIRwEVRwZSW4ohshAXe9LGG5SEAEQhQFCdmJ4ZoFjRIk4SK5kpCsnN8fp3pS\n0+nu6Znpnj4z+byep5+Zrj5V9a3TVfWtOlVdRyEEzMzMGm27RgdgZmYGTkhmZpYIJyQzM0tCj01I\nkj4uKeRef5P0vKSfStq/qOwMScleLMsty9BOjBskzah5UD2EpAlZHUxodCyNJmm+pF83Oo72SNpT\n0m2S/pp9d2d1YVozsmn0qWWMqZG0vaTPSHpY0mpJL0j6nqRdGh1bLfWGL/GjwHKgCdgXuBC4V9IB\nIYTVWZl/B+5qUHzV+DnwLuCFRgdi1g2mA+OBjxPX+WWNDKaHeDPwZeA7wELgAOCbwC7AyQ2Mq6Z6\nQ0JqCSE8k/3/oKS/APcA7wbmAYQQlhOTVpJCCC8BLzU6DrP2SNohhPB6FyczAlgcQrilFjFtI1YA\nI0II67P38yUNJyb1XqPHNtlVsCb7u31hQKkmu+w0/zJJn5f0R0lrJT0g6YCicpL0RUlPS9qUnSrP\nLD5Vzk3vHEl/kvSapJ9L2iN7/TQ71X5O0peLxt2qyU7SSZLuk/SSpHWS/kfS1M5WiqRhkn6UTe91\nSS2Sjs99vrOkpyT9TlK+7o6W9Iakz+SGfVbSb7Mml1WSHpL0waL5Dc2W6QxJV0hakdXxTZJ2krSf\npLuzZXumeNlyTTEHSbo/q88XJF0iqd31VtKHs7hey2L8maR9Olt/ZeaxLFuekyQ9KWm9pIWS/r6o\n3HxJ88uMf0PufWE9eHe2vqyVtFLS+dnnx2TrwXpJj0g6tExcx0l6LPuen5J0YokyoxSbzV6VtEHS\ng5LeU1TmBknLJb1L0m8kbSAelZerj4rbSmGdACYA79GW5vahFab5Jknfzbab17O/N0raoajosGx7\nW6e4/U3PryeS+kn6dlYv67L18XbFnXp+foXvYFy2vayR9BdJ35HUr6js2yTdma1jL0r6lqR/LrVM\n2fDFkjZKelnSLEm7FZX5QrYebci+l4XKttEQwsZcMirYH/hrubrrkUIIPfJFPDIIxC+lD7AD8cjr\nl8BKYJdc2RlxUduMH4hNBXcDk4GPAH8EngH65MpdnpWdCXwA+CKwDvgVsF3R9P5EbH77IPAJYnK8\nC3gQ+CpwJHBdVnZiiWUZmht2AfBp4OhsvEuAzcAZJZZjRjt1tTfwIvAYcGq2HD8E3gAm58odDLwO\nfD17P5h4ZHZr0fSuAj4JvD+b1swsjmNyZYbm6mR2ru42A3OAJcDngaOAuVksBxR/Z8CzwFeyevhW\n8fISd24BmJAbdkY27IfAROAfgSez77d/DdfBZdnyPZKtPx8C/gdYBQzIlZsPzC8z/g0l1oP/JTY9\n59eXb2R1dlI2nyeA54C+RfNZkcV0GnE9vCOr2yNy5Q4B1gO/zuKeCNyWffeH5srdAKzNpve5rK7/\nrkJ9VNxWiNvoOGAxsCj7fxywQ5npDczq4pVsWu8H/gn4ceF7zK0njwHnZHX2r9mw03LT2pXYdH8S\nsbnweGJLyqvAnmW+g0uy6V0I/A24OFeuL3HdXA5MzepwblZXxdvy14nr/beI6/FpwPPAw0BTVuYU\n4P+IzZlHZNObBnyyTN1Mz77XjzR6X1zLV8MD6HTgW1ac4tfzwNiisjMonZD+F9g+N+wj2fB3Z+93\nyzbSG4rGPTUrN7loektpm8z+JRv+1dywPsTkcH2JZRlaZlm3y8b7AbGpo3g5ZrRTV7OITYK7Fw2/\nh9jkmR/2xWzjO5KYTJ8HBlWYdiG2X5BLXGxJSPcVlZ+bDT81N2xgtjFeVPydAdOKxv8BcSc5IHs/\ngVxCApqB1cAPi8YbBmwCzqrhOriMuEMbmBs2Jovn5Nyw+XQsIU0vsb5sBoblhk/Oyo4vmk8AxuWG\nNQFPAb/KDbuXmKD7FpV7Evjv3LAbsukdV0VddGRb+XWp+igxzUuydfHgCmUK68lpRcOXAL+oMF4T\nsFO2Ln2xxHdwcVH5O4Cluff/nJU7LDdMxGTbui0Tt4O/5b/TbPjhWbl/yN7PBBZVud6dm4376Vqt\ny6m8ekOT3fHAWOAw4B+IR453ShpRxbj3hBA2594vyf4WmnbGEY+Ebioa78fEHej4EtP7v9z7p7K/\ndxcGZJ8/QzxrKUvS2yX9p6TniTujzcCniGeEHXUMcCewWlKfwiuLa5TaNj9enQ2/g3g0NyWE8HJR\nbIdKukPSSmI9bCae6ZSKbV7R+1J18ipxp1uqTn5a9P7HxKRzYJllfRfxQu+Pipb1uWze7y0zHpK2\ny48jqalc2ZzfZvEXFK9DndFaZ7n1ZWkI4Y+5MoV6LK6z50IID+XG/xvwM+CwbPl2JK63PwPeyNWP\niK0LxfWzmbgutKej20o1jgYeCSH8TxVlf170/jGKvgNJJyrepbYqi2k9cV0qtd4WT29J0fTGAX8O\nIfyuMCDEbHFz0XhHEQ/aitfHh4nJsFDfjwCjJV0j6UhJO5VaSEmDgSuAq0MI3y1VpifrDQnpsRDC\nwhDCIyGEW4lHjiIeObWnuP21cLG20FZcaONtc/dbtpN4Jfd5watF7zdVGN6PMiQ1E89eRhFP299D\nTLo/JDZ7dNQewBS2JLbC68rs890LBbON6sZsPotDCPcWxbY38Qh7N2Izzruz2O4qs0xdrZOVZd6/\nuURZiMsKcedavLwHkVvWEqYXlb+3QtmCNutQ2HLBv+z3W4VSdVOuHovnU1xfhWF9gTcRv7cmYjNU\ncf18FhiottfoXsqSWns6uq1UY3eqvxmp1LbcWjeSJgE/IZ4Fngz8HXG9fYnS31Wp6eW3vSHEg6hi\nxfVfWB+fYev67s+W9XEOcGYW193AXyXNLXF97e3Es+aU7xrutN5wl10bIYQNkv4AjKzB5Aor5Z7A\n44WB2RHO7tTvguK7gLcC7wkhtP6uRJ3/rcUrxHb8b5T5/C+5eexJbINfBBws6QshhH/NlT2G2B5/\nYoh3LxbGK3lEVwODgT8UvYfYlFjKK9nfj5P7znLWVpjX92l7NlCpbEdsJJ61FevMTro9g8sM20Tc\n+e5IvPZwLXEnuJUQwhv5t1XOtx7bysuUP/DoqJOAZ0IIHy8MULx5p7PfwQvAO0sML67/wvp4NFsf\nVLR+nh0IXgdcJ2kgW66Z/oSYpAo2AU+z5eatXqXXJaRsx7gvpXdGHfUQcQU4ibZHy/9IrLv5NZhH\nKYWde2tzYraSHtfJ6d1FTHKPhxA2lCskScQbEF5ny8Xcb0i6P4Tw+wqxvYPYJl6PW+tPJF4ULjiJ\neKF8Seni/IaYSPYLIczuyIxCCH8hl5xr6E/ACZL6hhA2AUh6L/EIudb2ljSu0GyXNTt+FPhdlmjW\nS/oV8ex7UVHy6Yp6bCu/AL4qaVQIYXEX49uJ2EyX9zHi2WJnPAScJumwQrNdtv2cUFTuHuIBwD4h\nhHuqmXDWBPwTSX8HnF702e+A4SVH7AV6Q0IaLWkQsZluCLHZYTfgmq5OOITwV0nfAs6XtJ54HWYE\ncBnxwmxxO3Ot/IZ4BHStpIuAnYl36b1MPDvpqOnA74AFkmYSL6YPJF6HeVsI4RNZubOJieh9IYRX\nJU0j3jTwn5LGZMnsl8QNe05WN0OAi4E/U58m4P+XNSE9Qrxz61PEmzhWlyocQlgj6UvEunsT8XrM\nauKR9njixfT/qEOclfyYeBH8h4q3eQ8j1nXJZeiilcSd2UXEM6IzgXdkfwvOBhYAd0uaRTzaH0S8\n+64phDCtozOt07bybWLz2i8lXUY8CBlEPDA7I4TQkTPYu4B/kPRt4lnwGGKT86pOxAXxho8vA3Ml\nfYVY158iblcQkxAhhGclfQOYqfgEmQeIZ8x7E68v/XsI4X5J3yceSP2W2BT4DmLC/EV+ppLGExP+\n0SGE+zoZe7J6Q0L6We7/l4gXM48JIdxdpnxHFVa2M4i3Yb9CbOo4v4ZHl22EEF7Kfn/wLeC/iEft\n/0pMtBd1Ynp/ljSGeF3tcuK1hFeIdTUbQNIh2WdXhBAeyMbbJOmfiM13/wKcGUJ4XNIpxDugbiPe\n+jqN2JQ3oZOLXMlxxIOLC4k78MuASyuNEEK4TtJzwJeIO7Q+xCa+XwEtdYixomyHcwbx7qgTiLeG\nn8rWF8Br4Rnib4UuJ15vWAb8Uwjh/lw8iySNJa5L3yEe5LxE/J6/14V513RbCSGsknQ48TufRmz6\nWwncx5ZraNX6ATEJfIJ41vEIMAno1I9zs23jaOK6+T3iWft/EG9W+Dq5g40QwgWSngQ+k70C8Sab\ne4l3+kL8achpxCS0K3Gbv4mtt3cRz+p6w/X/rSi7jdAsKYrP57uIeFt+cVOLWZIk3UF8osK+jY6l\nJ+oNZ0hmZt1O0tnEM6P/JV4P/Cjxx8hnVhrPynNCMjPrnNeJPyTfh9iM9jTwqRDCrIZG1YO5yc7M\nzJLQKy+MmZlZz+OEZGZmSXBCMjOzJDghNYC27n59bdZXymc78nggxT5a7sj+H140zXKvn5caP3WS\n9pI0W7EvmbWSfiJpQO7zsyQtURV9JXVy/snVb1fqJIv1452YZ13ruWhee0v6L8V+xNYoPtutqofW\nVjOupA8o9jm2QrGvpeWK/VCVeiSQdYdGP258W3yx5RH3HyE+Nfho4g/3AnBJldPYl/jjwDHZ+wFs\n6V9mHPGJFYHYr1J++NBS46f8Ij7ZYDnxR4wfIP6odBVwU67MjsS+gE6rUwxJ1W9X6ySL/eOdmG9d\n6zk3n52It1M/RnyK/3HEJzU8C+xci3GJfStdmW2H44k/Sn2c+JSUtzZ6vd8WXw0PYFt85RLSfkXD\n7wdWVxhvh9z/1xAfzV+u7OeyebyzzOcVx0/lRfxl+kPEx70oN/wSsic654Z9k/i8vq7Os18VZRpW\nv7Wok84mpFrWczvz+AKxH6H9csOGER9bdXYdx90/q5tz6rl8fpV+uckuLY8Auyh2eV7owvtAZV19\nk/UNpNh986nER5WUM5r4zKyniz8oN76kX0h6qET5gyRtzh4Z1N2OJz7t+OyQ7TEyfyZ2qbBXbtiP\ngXdKendHZqDYT9Bhki7Klr+a5641sn5rXicdjK1T9dxBk4GHQgjPFAaE2B/Ug7T/kOGujFt4Oref\nDtIATkhpGUY8sluXG3Yr8YGMk4kPm4TYNDSA+Gy2ckYR+4oq1ZdNufEfJHY50drviyQB3wV+E0L4\nUfWLUjOfID5w8g9q28FZc/Z5fsfRQnxA5THtTVTSIEknS7qR2AT1cDavFuBrVcTVyPqtR510JLaK\n01TUp4pXpSdtH0Bsciv2OKW7fej0uJKaJPWV9HZiFxArgP9sZx5WB35SQ2M1ZTuS/sRuFj4M3B5C\neC3uCwD4TmjbHxHEHV4Afk8J2TQPAMrt4MqN/yDxCPtgYpMQxI79xmXDupWkvsARxGsCm0sU2Uyu\nu4gQwhuSFhPjLTW9NxEf+nkssXO2vxGfRH0lcGcIoaouSxpZv7Wuk87EVsU0xxObn9vzAOUfyLsb\npfsP+itbnqhdTkfHfRg4NPv/GeLT7kt1vmd15oTUWE/l/n+DuIM7q6hMqacR7wWsCVnfOiXsT+wF\ns9yTrcuN/xBxJz0OeCi7Y+ubwMwQQqkjznp7J3HH+xli9xl5PwZeDVs/ePUl4qP7SxnFlp6EnyQ+\nQXpeaNuNfTUaWb+1rpPOxlZpmo8SE357atUBYld9jNiB4tuIT2S/R9LfhxCWNTSqbZATUmMdT7xT\nai3wpxDCxhJlXigxrB9bulsvZXT2t1ynZiXHDyGsKzry/RoxUZbs8kJSnxI7v1oamv39ddjSQSCS\nBhObN0s1q2wg3glWygPEOxqPzV63Amsk/ZLYb9JdIdcLbgXdUr9lDM3+1qpOOhtbpWmuo7puPio9\nt+xVSp/NlDv76fS4IYQns38fljSP2GXHNGI3GtaNfA2psR4LISwMITxdJhlB6Y32FeI1inJGZeOV\n22FWGv9BYJxi/0hnAF8KIbR2l5zdaHGRpEeIfTQhaYykeyUtlNQiaUpR+RmSFkl6RtLpxTOsoHDA\nVHyd5mPZ8t1QYpzdiB0ZbiWEsDmEcE8I4ewQwgjiEfEFwA7Zsjwn6feSPt9OXHWr3yrUtE66EFul\naY4nNh2297q3zPgQr/ccUGL4O4EnKi9G58cNIawiNtvt1848rA6ckHqmp4C+kt5S5vPRwLIKO5NK\n4/8aeCuxY7UHQwg3lZpACGFsCOEzWdPOLGBqCGEM8ZrAdMXeMQv6hhAOIe6oLpF0YDvLV7As+9u6\nc5G0J7Gnzu+HEJ4tMc4wStz5VmYZ/hhCuDaE8CFi52/HEDt/G9nOqHWv3wqWZX/rUScdia3SNAtN\ndu29Kh2c3EZMjm8rDJA0FDg8+6ySTo+bnWkOJ/5mybpbo+873xZflPkdUlGZGVmZPiU+G5p99uEy\n464E5laYdtnxiV19B+KdWiNLfB6At+TeTyT2jtmSey0Djs+VH5orfwPw+RKxzCgxLxGPaJ8l3mX4\nUeLO/rfATiXKDyA2M32qzHLvRNzZtPfap53vr2712111QonfIVUTWzX1XKNtZGfimcoS4q3ak4ln\npH8AmnPlxmfxTunEuLcQeyI+jnijyOlZXa4C3lGvZfOr/MvXkHqgEMIySb8jdsE8N/9ZdrS8B+Wb\nkyqOT2z/3wT8W8hdoyhRpnWWxB9JduQ3KflmyJ2zvytKxBkUu3L/AfE3WC8Tu3W+NITwWonpfjCL\nvVy31O8G7qkivrJ3f3VT/XZnnXQ0to5Os1NCCOslvY/4U4cbievZvcBZIYTi9a9Nl94dGPch4t2t\n5xDvMHwOmA9cEXxDQ2M0OiP61bkX8SxrNSWOirsyPvAt4o0Uu5YZLwADcu8HZuWPyQ17J9mRaFb+\nsuz/IcSd7IG5sv9MvGOrU8tRFNs84MZGfzddqd/uqhNKnyG1G1tq9exX73r5GlLPdRPx9yaf7ur4\nknaS9C5J5xEfu/LpEMLqaiYSQniVeMT8JcUHxD4BfId4xFmwWdIi4g9FLwptbyMeD3w7lD66r5qk\n0cD7gIu7Mp0a6kr9dluddDS2BOvZehH3GNuDSRoHHBJC+G5Xxifeen4r8DyxueLaGsYYgIEh3r1U\nN5KOyeaTzC/su6N+25l/2TrJvpfTiD8WrTq2FOvZeg8nJKur7kpI1jGFhBRCuKHRsZgVOCGZmVkS\nfA3JzMyS4IRkZmZJcEIyM7MkOCGZmVkSOvqkBt8B0RETJsS/8+c3Mgozs0ZT+0V8hmRmZolwQjIz\nsyQ4IZmZWRKckMzMLAlOSGZmlgT3h2Rm27zNmzezfPlyNm7c2OhQeqSmpiYGDBjAoEGD2G67zp/n\nOCGZ2TZv+fLl9O/fn6FDhyJVdYeyZUIIbN68mZUrV7J8+XL22WefTk/LTXZmts3buHEju+++u5NR\nJ0iib9++vPnNb2b9+vVdmpYTkpkZOBl1UVea6grcZGfblE/OXliX6c6aOqYu0zXblvgMyczMkuCE\nZGZmSXBCMjPr4c4//3yuvvrqqsoedthhPP7443WOqHOckMzMEjZ79myam5tpbm6mX79+NDU1tb4f\nMGAAK1euZM6cOZx++ulVTe/cc89l+vTpdY66c5yQzMwSNnXqVNatW8e6deu44IIL+NCHPtT6ftWq\nVcyZM4eJEyey4447VjW9yZMnc//997NixYo6R95xTkhmZj1ES0sLo0aNajNs3rx5jB8/vvX9unXr\naGpq4oUXXmgd9thjjzFkyBDWrl1Lv379OPTQQ7n77ru7Le5qOSGZmfUQLS0tjB49us2wJUuWsP/+\n+7e+b25uZvjw4SxatKh12LRp07jgggvo378/ACNGjGDx4sXdE3QH+HdIZmZ5Z50FLS31ncfo0VDl\nTQgFa9asYdmyZVslpFWrVrUmmoKxY8eyaNEiPvjBD7JgwQKeeOIJ5s6d2/p5//7925xBpcJnSGZm\nPcDixYvp378/w4YNazN84MCBrF27ts2wQkICOO+887j00kvp27dv6+dr165lwIAB9Q+6g3yGZGaW\n18Ezl+7S0tLCyJEjt3rE0ciRI1m6dCljx45tHTZ27Fi++c1vcvPNN7Nx40ZOPvnkNuM8+eSTnHrq\nqd0Sd0f4DMnMrAcodf0IYOLEiTzwwANtho0aNYoVK1ZwzjnncMUVV7RJYhs3buTRRx/lqKOOqnvM\nHeWEZGbWAyxevLhkQpoyZQp33nknGzZsaB22ww47cNBBBzF06FCOPfbYNuVvv/12JkyYwF577VX3\nmDvKTXZmZj3AwoWlHww8aNAgpkyZwnXXXcdZZ50FwKZNm3jxxReZOXPmVuWvuuoqZs2aVddYO8sJ\nycysh7v88svbvL/44os5/PDDGTdu3FZlH3744e4Kq8PcZGdm1kssWrSIXXfdlQULFnDNNdc0OpwO\n8xmSmVkvccghh7B69epGh9FpPkMyM7MkOCGZmVkSnJDMzCwJTkhmZpYEJyQzM0uCE5KZmSXBt32b\nJeyTs0v/Or8WZk0dU7dpm3WGz5DMzCwJPkMyq4F6nsmYtef8889n8ODBrc+yq4fDDjuM66+/ngMO\nOKBu8/AZkplZwmbPnk1zczPNzc3069ePpqam1vcDBgxg5cqVzJkzh9NPP72ucZx77rlMnz69rvPw\nGZKZWZF6n/F25Prd1KlTmTp1KgCXXHIJjz76KLfeemvr51deeSUTJ05kxx13rHmceZMnT+aMM85g\nxYoV7LnnnnWZh8+QzMx6iJaWFkaNGtVm2Lx58xg/fnzr+3Xr1tHU1MQLL7zQOuyxxx5jyJAhW3V1\n3hH9+vXj0EMP5e677+70NNrjhGRm1kOU6jV2yZIl7L///q3vm5ubGT58OIsWLWodNm3aNC644AL6\n9+/fpfmPGDGCxYsXd2kalTghmZn1AGvWrGHZsmVbJaRVq1ZtlWjGjh3bmpAWLFjAE088UZNrTP37\n92fVqlVdnk45TkhmZj3A4sWL6d+/P8OGDWszfODAgVs1xeUT0nnnncell15K375925R54403OhzD\n2rVrGTBgQIfHq5YTkplZD9DS0sLIkSOR1Gb4yJEjWbp0aZthhYR08803s3HjRk4++WQA5s+fz6RJ\nkzj++OO5/vrr+dznPscRRxzBkUceyfLly5k/fz5HH300kyZNYuzYsSxZsqTNdJ988smtrmHVkhOS\nmVkPUOr6EcDEiRN54IEH2gwbNWoUK1as4JxzzuGKK65ok8RWr17N3LlzGTx4MAMHDuT+++/na1/7\nGl//+tcBeO2117jtttuYM2cOX/nKV1rH27hxI48++ihHHXVUnZbQt32bmfUIixcv5swzz9xq+JQp\nUxg9ejQbNmxovfV7hx124KCDDqK5uZljjz22TfkxY8YgiSeeeIJbbrmFBQsWEEJg7733BuDggw9G\nEiNGjGhzp97tt9/OhAkT2Guvveq2jE5IZmY9wMKFpX8bNWjQIKZMmcJ1113X+qSGTZs28eKLLzJz\n5sytym+3XWwYGz58OCeeeCIXXnghAJs3b+bBBx+kpaWFEAJLly5lyJAhreNdddVVzJo1q9aL1YYT\nkplZkZ724NnLL7+8zfuLL76Yww8/nHHjxpUdZ9KkSdx3330cccQRSOKUU05h3333Zdddd2XSpEms\nXLmyTQJ6+OGH6xZ/gROSmVkvsWjRIo444ghGjhzJLbfcstXnEyZMYMKECQBI4uqrr27z+fz58xk+\nfDhXXXVVd4S7FSckM7Ne4pBDDmH16tWNDqPTnJDMzAxoewbVCL7t28zMkuCEZGZmSXBCMjOzJDgh\nmZlZEpyQzMyAEEKjQ+jRalF/Tkhmts1rampi8+bNjQ6jR9uwYQPbb799l6bhhGRm27wBAwawcuXK\nTnXJsK0LIfDaa6/x/PPPs8cee3RpWv4dkplt8wYNGsTy5ct5+umnGx1Kj7T99tszePBgdtllly5N\nxwnJzLZ52223Hfvss0+jw9jmucnOzMyS4IRkZmZJcEIyM7MkOCGZmVkSnJDMzCwJTkhmZpYEJyQz\nM0uCE5KZmSXBCcnMzJLghGRmZklwQjIzsyQ4IZmZWRKckMzMLAlOSGZmlgQnJDMzS4ITkpmZJcEd\n9FmXfHL2wrpMd9bUMXWZrpmly2dIZmaWBCckMzNLghOSmZklwQnJzMyS4IRkZmZJ8F12lqR63b1n\nZunyGZKZmSXBCcnMzJLghGRmZklwQjIzsyQ4IZmZWRKckMzMLAlOSGZmlgQnJDMzS4ITkpmZJcEJ\nyczMkuCEZGZmSXBCMjOzJDghmZlZEpyQzMwsCU5IZmaWBPeHZLaNqlefU7OmjqnLdK338xmSmZkl\nwQnJzMyS4IRkZmZJcEIyM7MkOCGZmVkSnJDMzCwJTkhmZpYEJyQzM0uCE5KZmSXBCcnMzJLghGRm\nZklwQjIzsyQ4IZmZWRKckMzMLAlOSGZmlgQnJDMzS4ITkpmZJcEJyczMkuCEZGZmSXBCMjOzJDgh\nmZlZEpyQzMwsCU5IZmaWBCckMzNLghOSmZklwQnJzMyS4IRkZmZJcEIyM7MkOCGZmVkSnJDMzCwJ\nTkhmZpYEJyQzM0uCE5KZmSXBCcnMzJLghGRmZklwQjIzsyT0aXQAZta7fHL2wrpMd9bUMXWZrqXD\nZ0hmZpYEJyQzM0uCE5KZmSXBCcnMzJLghGRmZklwQjIzsyQ4IZmZWRKckMzMLAlOSGZmlgQnJDMz\nS4ITkpmZJcEJyczMkuCEZGZmSXBCMjOzJDghmZlZEpyQzMwsCe6gbxtRr07TzMxqxWdIZmaWBCck\nMzNLgpvsEuJmNTPblvkMyczMkuCEZGZmSXBCMjOzJDghmZlZEpyQzMwsCU5IZmaWBCckMzNLghOS\nmZklwQnJzMyS4IRkZmZJcEIyM7MkOCGZmVkSnJDMzCwJTkhmZpYEJyQzM0uC+0Mysx6hXv2FzZo6\npi7TtY6wrw6+AAAGjUlEQVTzGZKZmSXBCcnMzJLghGRmZklwQjIzsyQ4IZmZWRKckMzMLAlOSGZm\nlgQnJDMzS4ITkpmZJcEJyczMkuCEZGZmSXBCMjOzJDghmZlZEpyQzMwsCU5IZmaWBCckMzNLghOS\nmZkloVf3GOseJs3Meo5enZDMzNrjA9d0uMnOzMyS4DOkTqj2iOpLK9YCcGWdjsDMzHoTnyGZmVkS\nnJDMzCwJTkhmZpYEJyQzM0uCE5KZmSXBCcnMzJLghGRmZklwQjIzsyQ4IZmZWRL8pAYzszrwM/I6\nrkMJyRVsZmb1ksQZUr0SnZlZb1PP/WWjTw58DcnMzJKgEEL1haW7gEFVFh8EvNyZoBrE8dZfT4vZ\n8dZfT4u5p8ULacT8cgjhmPYKdSghdYSkhSGEHnNxyPHWX0+L2fHWX0+LuafFCz0rZjfZmZlZEpyQ\nzMwsCfVMSN+v47TrwfHWX0+L2fHWX0+LuafFCz0o5rpdQzIzM+sIN9mZmVkSnJDMzCwJNU9IknaT\ndIuk9ZL+JOnkWs+jViTtIGlWFudaSS2Sjm10XNWQ9HZJGyXd1OhYqiHpJElPZuvFs5Le0+iYypE0\nVNKdkl6VtELSTElJPNUEQNJnJS2U9LqkG4o+e7+kpyS9Jul+SW9tUJhtlItZ0jhJ90j6q6SXJP1M\n0pAGhlqIq2wd58pMlxQkHdnN4ZWKpdI6sZOk70p6WdJqSQsaFGa76nGGdC2wCRgMnAL8m6QD6jCf\nWugDPAeMB3YFvgr8VNLQBsZUrWuBRxodRDUkHQV8AzgN6A+8F/hDQ4Oq7LvAi8AQYDRx/fh0QyNq\n6y/AZcAP8wMlDQLmAhcCuwELgZ90e3SllYwZGEi86D4UeCuwFri+WyMrrVy8AEjaF/go8EJ3BlVB\npXi/T1wfRmR/v9iNcXVITY/6JO0MnAAcGEJYB/xa0m3Ax4BptZxXLYQQ1gMzcoPukPRH4FBgWSNi\nqoakk4BVwG+A/RocTjUuBi4JITyUvX++kcFUYRgwM4SwEViRPaEkmYOqEMJcAEljgLfkPvow8HgI\n4WfZ5zOAlyUNDyE81e2B5pSLOYQwL19O0kzgge6NbmsV6rjgWuDLxIOXhisXr6ThwGTgLSGENdng\nR7s/wurU+gzpHcD/hRCW5oYtJqGNuRJJg4nL8HijYylH0i7AJcDZjY6lGpKagDHAmyQ9I2l51gS2\nY6Njq+Bq4KSsqePNwLHAXQ2OqRoHELc3oPWA61l6yPaXeS8Jb38Akj4KvB5CuLPRsVThMOBPwMVZ\nk90SSSc0Oqhyap2QmoE1RcNWE5tpkiZpe+BHwOxGH02241JgVghheaMDqdJgYHvgI8B7iE1gBxOb\nR1O1gLgTXwMsJzZ9/XdDI6pOM3F7y+sR2x+ApJHAdOBLjY6lHEn9gcuBLzQ6liq9BTiQuB7sBXwW\nmC1pREOjKqPWCWkdsEvRsF2I7cLJkrQdcCPx2tdnGxxOWZJGA0cC3250LB2wIft7TQjhhRDCy8C/\nABMbGFNZ2bpwF/FazM7EB1MOJF4DS12P3P4AJO0HzAO+EEL4VaPjqWAGcGMIYVmD46jWBmAzcFkI\nYVMI4QHgfuDoxoZVWq0T0lKgj6S354aNIuFTcEkCZhGP5E8IIWxucEiVTCBe/P2zpBXAucAJkhY1\nMqhKQgivEs8y8r/ATvnX2LsB+xCvIb0eQniFeJE9yQRa5HHi9ga0XtPdl4S3P4DsTsBfApeGEG5s\ndDzteD/w+ezuyxXA3sQbob7c4LjK+X2JYclufzVNSFmb9VzgEkk7SzocOI549pGqfyPefTIphLCh\nvcIN9n3iDmZ09voe8HPgA40MqgrXA5+TtIekgcS7fO5ocEwlZWdwfwTOlNRH0gBgKqU37IbI4uoH\nNAFNkvplt6XfAhwo6YTs8+nA71Nogi4Xc3aN7j7iAcD3GhvlFhXq+P3EJrDCNvgX4HTiTQ4NUyHe\nBcCfgfOzMocDRwB3NzDc8kIINX0RjzD/G1hPrIiTaz2PGsb6VuLRwkZic0fhdUqjY6sy/hnATY2O\no4o4tyfejbQKWAF8B+jX6LgqxDsamA+8SuxH5qfA4EbHVfS9h6LXjOyzI4GniE0184GhjY63UszA\nRdn/+e1vXarxlii3DDgy5XiJ10N/m+2TnwCOb3S85V5+lp2ZmSXBjw4yM7MkOCGZmVkSnJDMzCwJ\nTkhmZpYEJyQzM0uCE5KZmSXBCcnMzJLghGRmZklwQjIzsyT8fxGbSa78glC/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f79a17b67b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "plt.hist(\n",
    "    Tyr,\n",
    "    np.arange(19),\n",
    "    align='left',\n",
    "    label='$T(y_\\mathrm{rep})$',\n",
    "    color=plot_tools.lighten('C0', 0.3)\n",
    ")\n",
    "plt.axvline(Ty, color='red', label='$T(y)$')\n",
    "plt.xlim((-0.5, 17.5))\n",
    "plt.title(\n",
    "    'Binomial example - number of changes? \\n'\n",
    "    r'$\\operatorname{Pr}(T(y_\\mathrm{rep},\\theta) \\leq T(y,\\theta)|y) = 0.03$',\n",
    "    fontsize=16\n",
    ")\n",
    "plt.legend()\n",
    "plot_tools.modify_axes.only_x(plt.gca())\n",
    "plt.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
